The Unique Features and Performance of DifferentialEquations.jl


JuliaCon 2017


Chris Rackauckas

University of California, Irvine

Quick About Me

  • PhD Candidate in Mathematics at UC Irvine
  • Research in:
    • Numerical methods for stochastic ODEs/PDEs
    • Multiscale and spatial biological simulations
    • Effects of stochasticity in gene regulatory networks
  • I created DiffEq because it was necessary:
    • No software was using the latest methods
    • I wanted the latest and fastest methods to benchmark against
    • I am too lazy to change my code around: everything needs one interface
    • I wanted to be able to use these methods to build large biological models
    • I wanted parallel computing to be part of the design
    • I wanted tools to be able to match models to data

DifferentialEquations.jl is about modeling

Scientific laws, theories, and experiments

tell us how things change


Differential equations give a mechanistic model


Differential equations map "change" to simulations

</center> Differential equations bridge from models to predictions

Differential Equations vs Machine Learning

  • The models are different
    • Machine learning uses general non-mechanistic models
    • Differential equations utilize knowledge about the structure of the problem
  • The reliance on data is different
    • Machine learning needs to learn from data
    • Machine learning models stack simple models together
    • DE models are complex and nonlinear
    • DE can model without data, or can calibrate to data
    • DE have fewer parameters, and can be learned with less data
  • Both are used to create predictions
    • Ordinary and partial differential equations are the standard tools of physics
    • Stochastic differential equations are a standard tool in stock market analysis

There are many types of differential equations

  • Ordinary differential equations evolve over one variable (let's call this time)
  • Systems of differential equations evolve multiple quantities at the same time
  • Stochastic differential equations add randomness (referred to as "noise")
  • Delay differential equations allow for delayed interactions
  • Gillespie simulations evolve discrete quantities via a stochastic model
  • Partial differential equations evolve over multiple variables (space and time)

DifferentialEquations.jl incorporates all of these under a single interface

Problem 1: Ordinary Differential Equations

http://docs.juliadiffeq.org/latest/tutorials/ode_example.html

$$\frac{du(t)}{dt} = f(t,u(t))$$

where $u(t_0) = u_0$ is the known initial value. We want to know $u(t)$ for $t\in[t_0,t_f]$

Exponential Growth Model

  • My money has been temporarily stolen by the banking cartel (it's in the bank)
  • My money grows at a rate of 98% per year due to interest (really good bank)
  • I start with $1: How much will I have in 5 years?
$$ \frac{du}{dt} = 0.98u $$$$u(0)=1$$

Solve for $u(t)$ for $t\in[0,5]$

In [1]:
using DifferentialEquations
f(t,u) = 0.98*u
u0=1.0
tspan = (0.0,5.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob);
In [2]:
using Plots; plotly()
# Use the plot recipe, along with some plotting attributes
plot(sol,linewidth=5,title="My Bank Account",
     xaxis="Time (t)",yaxis="Money!",legend=false)
Out[2]:

Controlling the Solvers

See http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html

In [3]:
sol = solve(prob)
println(length(sol.t))
sol = solve(prob,abstol=1e-6,reltol=1e-6)
println(length(sol.t))
sol = solve(prob,saveat=0.5,reltol=1e-6)
println(sol.t)
sol = solve(prob,tstops=[0.5],reltol=1e-4)
println(sol.t)
sol = solve(prob,reltol=1e-6,save_everystep=false)
println(sol.t)

11
15
[0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0]
[0.0,0.0634774,0.221501,0.43294,0.5,0.696266,0.923443,1.20696,1.53056,1.89898,2.30449,2.74501,3.21551,3.71248,4.23208,4.77109,5.0]

Analyzing the Solution

http://docs.juliadiffeq.org/latest/basics/solution.html

In [4]:
sol = solve(prob,Tsit5())
println(sol.t[5]) # 5th timepoint
println(sol[5]) # 5th timepoint value
println([t+2u for (t,u) in zip(sol.t,sol.u)])

1.106052911995544
2.956279956560924
[2.0,2.30727,3.17657,4.63952,7.01861,11.2071,18.9435,34.2493,66.3366,137.407,273.572]

DifferentialEquations.jl comes with interpolations on the solution type

In [5]:
sol(0.45) # The value of the solution at t=0.45
Out[5]:
1.5542610480525965

Choosing Solvers

http://docs.juliadiffeq.org/latest/solvers/ode_solve.html

DifferentialEquations.jl is composed of many (30+?) packages which together provide its functionality under one interface. For Ordinary Differential Equations, we have:

  • OrdinaryDiffEq.jl: native Julia library. Fully featured, and in many cases faster than the classic C/Fortran libraries.
  • Sundials.jl: provides Adams and BDF methods. A standard library
  • LSODA.jl: provides the stability-detecting LSODA algorithm
  • ODEInterface.jl: provides classic Fortran algorithms like dopri and radau
  • ODE.jl: backwards compatibility

DifferentialEquations.jl also provides automatic detection of the appropriate algorithm

In [6]:
# Automatically uses CVODE_BDF from Sundials
sol = solve(prob,alg_hints=[:stiff],reltol=1e-8,abstol=1e-8) 
# Use Tsit5 from OrdinaryDiffEq
sol = solve(prob,Tsit5()); 
# Use CVODE_BDF from Sundials
sol = solve(prob,CVODE_BDF()); 
using ODEInterfaceDiffEq
# Use radau from ODEInterfaceDiffEq
sol = solve(prob,radau());

The *DiffEq Libraries

The *DiffEq libraries (OrdinaryDiffEq.jl, StochasticDiffEq.jl, DelayDiffEq.jl) are the native Julia solvers. They have a lot of unique features:

  • They routinely benchmark as the fastest implementations
  • They have cheap/free high order interpolations (will show in a second!)
  • They have expressive event handling
  • They can handle generic Number and AbstractArray types

The wrapped libraries, while classics (Sundials.jl, ODEInterfaceDiffEq.jl) are restricted to Vector{Float64}.

http://docs.juliadiffeq.org/latest/basics/compatibility_chart.html

Method Recommendations

  • BS3() for fast low accuracy non-stiff
  • Tsit5() for standard non-stiff. This is the first algorithm to try in most cases.
  • Vern7() for high accuracy non-stiff
  • Rosenbrock23() for stiff equations with Julia-defined types, events, etc.
  • radau() for stiff equations at higher accuracy on Vector{Float64}
  • CVODE_BDF() for large stiff equations (discretizations of PDEs) on Vector{Float64}

Translations from MATLAB/Python/R

  • ode23 --> BS3()
  • ode45/dopri5 --> DP5(), though in most cases Tsit5() is more efficient
  • ode23s --> Rosenbrock23()
  • ode113 --> CVODE_Adams(), though in many cases Vern7() or Vern8() are more efficient
  • dop853 --> DP8(), though in most cases Vern7() or Vern8() are more efficient
  • ode15s/vode --> CVODE_BDF(), though in many cases radau() is more efficient
  • ode23t --> Trapezoid()
  • lsoda --> lsoda() (requires Pkg.add("LSODA"); using LSODA)
  • ode15i --> IDA()

Example 2: Systems of Differential Equations

$$ \begin{align} \frac{dx}{dt} &= σ(y-x) \\ \frac{dy}{dt} &= x(ρ-z) - y \\ \frac{dz}{dt} &= xy - βz \\ \end{align} $$
In [7]:
function lorenz(t,u,du)
 du[1] = 10.0(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob,Tsit5());

The Plot Recipe Interpolates and Does Phase Plots

In [8]:
plot(sol,vars=(1,2,3))
Out[8]:
In [9]:
plot(sol,vars=(1,2,3),denseplot=false)
Out[9]:

ParameterizedFunctions

$$ \begin{align} \frac{dx}{dt} &= σ(y-x) \\ \frac{dy}{dt} &= x(ρ-z) - y \\ \frac{dz}{dt} &= xy - βz \\ \end{align} $$
In [10]:
g = @ode_def_nohes LorenzExample begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ=>10.0 ρ=>28.0 β=(8/3)
Out[10]:
(::LorenzExample) (generic function with 15 methods)

DifferentialEquations.jl respects the user's input type

In [11]:
using StaticArrays, DifferentialEquations
A  = big.(@SMatrix [ 1.0  0.0 0.0 -5.0
                4.0 -2.0 4.0 -3.0
               -4.0  0.0 0.0  1.0
                5.0 -2.0 2.0  3.0])
u0 = big.(@SMatrix rand(4,2))
tspan = (0.0,1.0)
f_SA(t,u) = A*u
prob = ODEProblem(f_SA,u0,tspan)
sol = solve(prob,Tsit5())
using Plots; plotly(); plot(sol)
Out[11]:
In [12]:
using Unitful
f_units(t,y) = 0.5*y / 3.0u"s"
u0 = 1.0u"N"
prob = ODEProblem(f_units,u0,(0.0u"s",1.0u"s"))
sol =solve(prob,Tsit5());
println(sol.t[end])
println(sol[end])
sol(0.5u"s")

1.0 s
1.1813604129914028 N
In [13]:
using DifferentialEquations, SymEngine
y0 = symbols(:y0)
u0 = y0
f_sym(t,u) = 2u
prob = ODEProblem(f_sym,u0,(0.0,1.0))
sol = solve(prob,RK4(),dt=1/10);

sol = solve(prob,RK4(),dt=1/1000)
end_solution  = lambdify(sol[end])

println(end_solution(2)) # 14.778112197857341
println(2exp(2)) # 14.7781121978613
println(end_solution(3)) # 22.167168296786013
println(3exp(2)) # 22.16716829679195

14.

DifferentialEquations.jl is unique when speed matters

DifferentialEquations.jl is a research lab which includes many new methods

In [14]:
using DifferentialEquations, Plots, ODEInterfaceDiffEq, ODE
# 2D Linear ODE
function test_f(t,u,du)
  for i in 1:length(u)
    du[i] = 1.01*u[i]
  end
end
function test_f(::Type{Val{:analytic}},t,u₀)
  u₀*exp(1.01*t)
end
tspan = (0.0,10.0)
prob = ODEProblem(test_f,rand(100,100),tspan)

abstols = 1./10.^(3:13)
reltols = 1./10.^(0:10);
setups = [Dict(:alg=>DP5())
          Dict(:alg=>ode45())
          Dict(:alg=>dopri5())
          Dict(:alg=>Tsit5())]
names = ["DifferentialEquations";"ODE";"ODEInterface";"DifferentialEquations Tsit5"]
wp = ode_workprecision_set(prob,abstols,reltols,setups;names=names,save_everystep=false);
In [15]:
plot(wp)
Out[15]:

Callbacks and Event Handling

You can directly apply dynamic behavior to an ODE through its callbacks and event handling

In [16]:
function b_condition(t,u,integrator) # Event when event_f(t,u) == 0
  u[1]
end

function b_affect!(integrator)
  integrator.u[2] = -integrator.u[2]
end

cb = ContinuousCallback(b_condition,b_affect!)

bb = @ode_def_bare BallBounce begin
  dy =  v
  dv = -g
end g=9.81

u0 = [50.0,0.0]
tspan = (0.0,15.0)
prob = ODEProblem(bb,u0,tspan)
sol = solve(prob,Tsit5(),callback=cb);
In [17]:
plot(sol,plotdensity=1000)
Out[17]:

Pretty much anything can be changed with the integrator

Let's model a growing cell population

  • Each cell has an internal protein X
  • When the concentration of X hits 1, the cell divides
  • When there are more cells, there is less food and so X is created less rapidly
In [18]:
using DifferentialEquations, Plots
function growth(t,u,du)
  for i in 1:length(u)
    du[i] = 0.3u[i]/length(u)
  end
end

function condition(t,u,integrator) # Event when event_f(t,u) == 0
  1-maximum(u)
end

function affect!(integrator)
  u = integrator.u
  resize!(integrator,length(u)+1)
  maxidx = findmax(u)[2]
  Θ = rand()
  u[maxidx] = Θ
  u[end] = 1-Θ
  length(u) == 10 && terminate!(integrator)
  nothing
end

callback = ContinuousCallback(condition,affect!)
u0 = [0.2]
tspan = (0.0,100.0)
prob = ODEProblem(growth,u0,tspan)
sol = solve(prob,callback=callback);
In [19]:
plot(sol.t,map((x)->length(x),sol[:]),lw=3,
     ylabel="Number of Cells",xlabel="Time")
Out[19]:
In [20]:
ts = linspace(0,sol.t[end],1000)
plot(ts,map((x)->x[1],sol.(ts)),lw=3,
     ylabel="Amount of X in Cell 1",xlabel="Time")
Out[20]:

This flexibility takes DiffEq from equation solver to simluation engine

  • Abstract types can be used to encode complex models
  • Integrators and event handling gives dynamic interactions
  • Speed, control, adaptivity, plotting, etc. then all comes for free!

This leads to more sustainable scientific software development:

  • Many scientific projects develop their own internal solvers for their models
    • This requires a lot of extra time!
    • How well tested is it? How well optimized is it? How maintainable is it?
    • Is it using all of the latest PI-stepsize adaptivity? Verner's coefficients? etc.?
  • DifferentialEquations.jl is thus ideal as a core engine in complex models
    • Developer-friendly dependency management, MIT licensed
    • Active developer channels (Gitter, or anywhere...)

Example: Multiscale Biological Models

In [21]:
using MultiScaleArrays
using OrdinaryDiffEq, DiffEqBase, StochasticDiffEq

# Define a hierarchy

immutable Cell{B} <: AbstractMultiScaleArrayLeaf{B}
  x::Vector{B}
end
immutable Population{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
  x::Vector{T}
  y::Vector{B}
  end_idxs::Vector{Int}
end
immutable Tissue{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
  x::Vector{T}
  y::Vector{B}
  end_idxs::Vector{Int}
end
immutable Embryo{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
  x::Vector{T}
  y::Vector{B}
  end_idxs::Vector{Int}
end
In [22]:
cell1 = Cell([1.0;2.0;3.0])
cell2 = Cell([4.0;5])
population = construct(Population,deepcopy([cell1,cell2])) # Make a Population from cells
cell3 = Cell([3.0;2.0;5.0])
cell4 = Cell([4.0;6])
population2 = construct(Population,deepcopy([cell3,cell4]))
tissue1 = construct(Tissue,deepcopy([population,population2])) # Make a Tissue from Populations
tissue2 = construct(Tissue,deepcopy([population2,population]))
embryo = construct(Embryo,deepcopy([tissue1,tissue2])); # Make an embryo from Tissues
In [23]:
function cell_ode(t,cell,dcell)
  m = mean(cell)
  for i in eachindex(cell)
    dcell[i] = -m*cell[i]
  end
end
function embryo_ode(t,embryo,dembryo)
  for (cell,y,z) in LevelIterIdx(embryo,2)
    cell_ode(t,cell,@view dembryo[y:z])
  end
end

# Have a cell divide at t=0.5
tstop = [0.5]
function grow_condition(t,u,integrator)
  t ∈ tstop
end
function grow_affect!(integrator)
  add_daughter!(integrator,integrator.u[1,1,1],1,1)
end
growing_cb = DiscreteCallback(grow_condition,grow_affect!)

prob = ODEProblem(embryo_ode,embryo,(0.0,1.0));
sol = solve(prob,Tsit5(),callback=growing_cb,tstops=tstop)

sol[end]
Out[23]:
Embryo{Tissue{Population{Cell{Float64},Float64},Float64},Float64}(Tissue{Population{Cell{Float64},Float64},Float64}[Tissue{Population{Cell{Float64},Float64},Float64}(Population{Cell{Float64},Float64}[Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.290911,0.581821,0.872732]),Cell{Float64}([1.16364,1.45455]),Cell{Float64}([0.0,0.0,0.0])],Float64[],[3,5,8]),Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.600023,0.400015,1.00004]),Cell{Float64}([0.80003,1.20005])],Float64[],[3,5])],Float64[],[8,13]),Tissue{Population{Cell{Float64},Float64},Float64}(Population{Cell{Float64},Float64}[Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.600023,0.400015,1.00004]),Cell{Float64}([0.80003,1.20005])],Float64[],[3,5]),Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.250004,0.500008,0.750013]),Cell{Float64}([1.00002,1.25002])],Float64[],[3,5])],Float64[],[5,10])],Float64[],[13,23])

That's the simulation of deterministic models

Let's add stochasticity

Stochastic Differential Equations

$$ du = f(t,u)dt + \sum_i g_i(t,u)dW^i $$
  • The $f$ term is the deterministic or "drift" term
  • The $g$ term is the noise or "diffusion term

Well-known model: the Black-Scholes Model

  • My stock is growing at an average rate of 101%
  • However, the value of my money has fluctuatations of 87% of its current value
$$ du = 1.01u + 0.87dW_t $$
In [24]:
using DifferentialEquations, Plots; plotly()
α=1.01; β=0.87; u₀=1.0
bs_f(t,u) = α*u
bs_g(t,u) = β*u
bs_f(::Type{Val{:analytic}},t,u₀,W) = u₀*exp((α-(β^2)/2)*t+β*W)
prob = SDEProblem(bs_f,bs_g,u₀,(0.0,1.0))
sol = solve(prob,EM(),dt=1/2^(4))
plot(sol,plot_analytic=true)
Out[24]:

Fast adaptive methods are unique to DifferentialEquations.jl

In [25]:
sol =solve(prob,SRIW1())
plot(sol,plot_analytic=true)
Out[25]:

They tend to be orders of magnitude faster than standard methods

Many models cannot be solved without these new methods. Even those that can are faster with them.

Gillespie Simulations

Gillespie simulations are discrete models which change stochastically over time.

  • A protein binds to DNA at a rate X
  • A protein is created at a rate Y

These types of models are related to differential equations

  • When the number of objects gets high, it converges for an SDE on the concentrations
  • When the number gets really high, it converges to a deterministic ODE
In [38]:
prob = DiscreteProblem([999,1,0],(0.0,250.0))

rate = (t,u) -> (0.1/1000.0)*u[1]*u[2]
jump_affect! = function (integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
jump = ConstantRateJump(rate,jump_affect!)

rate = (t,u) -> 0.01u[2]
jump_affect! = function (integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,jump_affect!)
jump_prob = JumpProblem(prob,Direct(),jump,jump2)

sol = solve(jump_prob,Discrete(apply_map=false));
In [39]:
plot(sol)
Out[39]:

Mix jumps and differential equations

In [40]:
diff_f = function (t,u,du)
  du[4] = u[2]*u[3]/100000 - u[1]*u[2]/100000
end
diff_g = function (t,u,du)
  du[4] = 0.1u[4]
end

prob = SDEProblem(diff_f,diff_g,[999.0,1.0,0.0,1.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,SRIW1()); plot(sol)
Out[40]:

These features together allow for expressive scientific modeling

Let's go back to the multiscale embyro model:

  • We can have cell proteins change according to a stochastic differential equation
  • A rate of division can be dependent on a specific "growth factor"
  • "Jumps" using this rate can cause cell division, growing our multiscale model
  • We get high performance algorithms and optimized implementations the whole way down.

<3 Julia

Addons

Up until now, I have shown how to solve many different equations and encode a wide variety of models to simulate using the various solvers in DifferentialEquations.jl. But because all of this uses a single interface, an addon suite was able to be developed.

  • Tools for automatic parallelization of Monte Carlo simulations
  • Uncertainty quantification
  • Parameter estimation and data fitting

Monte Carlo Simulations

Many times you need to solve the same equation multiple times. This is especially true for stochastic simulations.

In [29]:
pf_func = function (t,u,p,du)
  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
  du[2] = -3 * u[2] + u[1]*u[2]
end
pf = ParameterizedFunction(pf_func,[1.5,1.0])
pg_func = function (t,u,p,du)
  du[1] = p[1]*u[1]
  du[2] = p[2]*u[2]
end
pg = ParameterizedFunction(pg_func,[0.1,0.1])
prob = SDEProblem(pf,pg,[1.0,1.0],(0.0,10.0))

prob_func = function (prob,i)
  set_param_values!(prob.g,0.3rand(2))
  prob
end
monte_prob = MonteCarloProblem(prob,prob_func=prob_func)
sim = solve(monte_prob,SRIW1(),num_monte=10);
In [30]:
plot(sim,linealpha=0.6,color=:blue,vars=(0,1),title="The Timeseries'")
plot!(sim,linealpha=0.6,color=:red,vars=(0,2))
Out[30]:
In [31]:
summ = MonteCarloSummary(sim,0:0.1:10)
gr() # Note that plotly does not support ribbon plots
plot(summ,fillalpha=0.4)
Out[31]:
0 2 4 6 8 10 0 5 10

Uncertainty Quantification

Estimate the uncertainty in the solution due to numerical error

In [32]:
fitz = @ode_def_nohes FitzhughNagumo begin
  dV = c*(V - V^3/3 + R)
  dR = -(1/c)*(V -  a - b*R)
end a=0.2 b=0.2 c=3.0
u0 = [-1.0;1.0]
tspan = (0.0,20.0)
prob = ODEProblem(fitz,u0,tspan)
cb = AdaptiveProbIntsUncertainty(5)
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Tsit5(),num_monte=100,callback=cb,abstol=1e-3,reltol=1e-1)
using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)
Out[32]:
In [33]:
sim = solve(monte_prob,Tsit5(),num_monte=100,callback=cb,abstol=1e-6,reltol=1e-3)
plot(sim,vars=(0,1),linealpha=0.4)
Out[33]:

Case Study: Parameter Estimation

Instead of simply solving differential equation models, one can ask the inverse problem: given that I see data like y, what needs to be true about the model such that we can match the data?

The act of estimating the parameters for a mechanistic differential equation model is called parameter estimation.

  • Parameter estimation required repeated solving of differential equations
    • Solving equations must be fast
  • It requires access to optimization routines (global optimization is necessary)
  • Machine learning tools (regularization, cross validation, etc.) are a must

Thus DifferentialEquations.jl is well-suited for handling these kinds of problems.

Problem Setup

In [34]:
using DifferentialEquations, Plots; plotly()
lotka = @ode_def_nohes LotkaVolterraTest begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a=>1.5 b=>1.0 c=>3.0 d=>1.0

u0 = [1.0;1.0]; tspan = (0.0,10.0); prob = ODEProblem(lotka,u0,tspan)
sol = solve(prob,Vern8(),abstol=1e-14,reltol=1e-14); t = collect(linspace(0,10,200))
randomized = [(sol(t[i]) + .01randn(2)) for i in 1:length(t)]
using RecursiveArrayTools; data = vecvec_to_mat(randomized); 
plot(sol,plotdensity=1000); scatter!(t,data)
Out[34]:

We can tell it to build an objective function that internally uses the DiffEq solver

This objective function has dispatches to work with all of the common optimization packages in Julia (JuMP/MathProgBase, Optim, BlackBoxOptim, ...)

In [35]:
obj = build_loss_objective(prob,Tsit5(),CostVData(t,data),maxiters=10000,verbose=true)
Out[35]:
(::DiffEqObjective) (generic function with 2 methods)

Let's recover the parameters using the NLopt.jl package

In [36]:
using NLopt
opt = Opt(:LN_BOBYQA, 4)
lower_bounds!(opt,zeros(4))
upper_bounds!(opt,5ones(4))
min_objective!(opt, obj)
maxeval!(opt, Int(1e5))
(minf,minx,ret) = NLopt.optimize(opt,ones(4))

WARNING: using NLopt.population in module Main conflicts with an existing identifier.
Out[36]:
(0.009874582531472897,[1.51121,1.00488,2.9638,0.987147],:ROUNDOFF_LIMITED)

What's Next?

The next focus for DifferentialEquations.jl development is PDEs and stiff problems

  • Native Julia stiff solvers are currently a lower order of accuracy then the top methods (@miguelraz's GSoC project will help here)
  • Parallel-in-time Neural Network methods for large systems of stiff ODEs/SDEs (@akaysh's GSoC project)
  • Fast matrix-free operators for FDM derivative discretizations (@shivin9's GSoC project)
  • A wrapper for FEniCS's FEM toolbox (@ysimillides's GSoC project)
  • High order Rosenbrock methods
  • Stiffness detection and switching (already implemented, but no default methods)

Other areas which are under development are:

  • New high order (+ adaptive) methods for SDEs will be published soon (orders of magnitude improvement)
  • Integration of regularization, maximum likelihood, and Bayesian (MCMC) techniques for parameter estimation (@ayush-iitkgp's GSoC project)
  • Boundary value problem solvers (@YingboMa's GSoC project)
  • More geometric and symplectic integrators
  • Integration of DiffEq with applications (large scale biological models, Pk/Pd estimation, tools for physical and financial modeling)